Hasta ahora hemos encontrado clusters indicando cuantos se necesitaban, e indirectamente hemos forzado a que cada caso sea parte de uno de esos clusters. Veamos la data nuevamente:
# bibliotecas:
library(stringr)
library(magrittr)
library(htmltab)
library(factoextra)
library(cluster)
# coleccion
links=list(web="https://en.wikipedia.org/wiki/Democracy_Index",
xpath ='//*[@id="mw-content-text"]/div/table[2]/tbody')
demo<- htmltab(doc = links$web, which =links$xpath)
# limpieza
names(demo)=str_split(names(demo),">>",simplify = T)[,1]%>%gsub('\\s','',.)
demo[,-c(1,8,9)]=lapply(demo[,-c(1,8,9)], trimws,whitespace = "[\\h\\v]")
# preparación
demo=demo[,-c(1)] #sin Rank
demo[,-c(1,8,9)]=lapply(demo[,-c(1,8,9)], as.numeric) # a numerico
row.names(demo)=demo$Country # cambiando row.names
demo=na.omit(demo)
# veamos que tenemos:
str(demo)
## 'data.frame': 167 obs. of 9 variables:
## $ Country : chr "Norway" "Iceland" "Sweden" "New Zealand" ...
## $ Score : num 9.87 9.58 9.39 9.26 9.22 9.15 9.15 9.14 9.09 9.03 ...
## $ Electoralprocessandpluralism: num 10 10 9.58 10 10 9.58 9.58 10 10 9.58 ...
## $ Functioningofgovernment : num 9.64 9.29 9.64 9.29 9.29 7.86 9.64 8.93 8.93 9.29 ...
## $ Politicalparticipation : num 10 8.89 8.33 8.89 8.33 8.33 7.78 8.33 7.78 7.78 ...
## $ Politicalculture : num 10 10 10 8.13 9.38 10 8.75 8.75 8.75 9.38 ...
## $ Civilliberties : num 9.71 9.71 9.41 10 9.12 10 10 9.71 10 9.12 ...
## $ Regimetype : chr "Full democracy" "Full democracy" "Full democracy" "Full democracy" ...
## $ Continent : chr "Europe" "Europe" "Europe" "Oceania" ...
Nuestro punto de partida clave siempre ha sido el cálculo de la matriz de distancias, añadamos la semilla aleatoria:
set.seed(2019)
inputData=demo[,c(3:7)]
g.dist = daisy(inputData, metric="gower")
Y con esa matriz calculamos cuatro clusters en nuestras clases previas, pero tal cantidad de clusters salió de una decisión algo arbitraria. Una pregunta exploratoria clave es cuántos clusters deberíamos calcular, y según ellos saber que hay una cantidad diferenciada de perfiles.
Existe la medida gap, que sirve para determinar el mejor numero de clusters a pedir. Veamos:
fviz_nbclust(inputData, pam,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)
fviz_nbclust(inputData, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)
El numero de clusters varía, pero quedemonos con seis en ambos enfoques.
Clusterizemos:
res.pam = pam(g.dist,6,cluster.only = F)
res.agnes = hcut(g.dist, k = 6,hc_func='agnes',hc_method = "ward.D")
res.diana = hcut(g.dist, k = 6,hc_func='diana')
Para evaluar, podemos analizar las siluetas (silhouettes), una medida que indica la calidad de asignación de un caso particular.
fviz_silhouette(res.pam)
## cluster size ave.sil.width
## 1 1 14 0.56
## 2 2 29 0.25
## 3 3 31 0.34
## 4 4 35 0.09
## 5 5 32 0.26
## 6 6 26 0.28
fviz_silhouette(res.agnes)
## cluster size ave.sil.width
## 1 1 14 0.57
## 2 2 33 0.26
## 3 3 34 0.25
## 4 4 34 0.08
## 5 5 19 0.34
## 6 6 33 0.25
fviz_silhouette(res.diana)
## cluster size ave.sil.width
## 1 1 39 0.47
## 2 2 52 0.27
## 3 3 3 0.25
## 4 4 38 0.10
## 5 5 18 0.34
## 6 6 17 0.25
Se asume que el gráfico que tiene menos siluetas negativas es el preferible a los demás.
Esta etapa es para identificar a los casos mal asignados: los que tienen silueta negativa. Para ello es bueno saber lo que cada resultado nos trajo:
# por ejemplo tiene:
str(res.pam)
## List of 9
## $ medoids : chr [1:6] "Finland" "Cape Verde" "Guyana" "Benin" ...
## $ id.med : int [1:6] 8 26 55 82 126 155
## $ clustering: Named int [1:167] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:167] "Norway" "Iceland" "Sweden" "New Zealand" ...
## $ objective : Named num [1:2] 0.0831 0.0794
## ..- attr(*, "names")= chr [1:2] "build" "swap"
## $ isolation : Factor w/ 3 levels "no","L","L*": 1 1 1 1 1 1
## ..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
## $ clusinfo : num [1:6, 1:5] 14 29 31 35 32 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:5] "size" "max_diss" "av_diss" "diameter" ...
## $ silinfo :List of 3
## ..$ widths : num [1:167, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:167] "Finland" "Iceland" "Denmark" "Sweden" ...
## .. .. ..$ : chr [1:3] "cluster" "neighbor" "sil_width"
## ..$ clus.avg.widths: num [1:6] 0.5571 0.2517 0.3405 0.0867 0.2638 ...
## ..$ avg.width : num 0.266
## $ diss : NULL
## $ call : language pam(x = g.dist, k = 6, cluster.only = F)
## - attr(*, "class")= chr [1:2] "pam" "partition"
Aquí sabemos que res.pam es una lista con varios elementos. Uno de ellos es la información de siluetas, el cual tiene otros componentes:
str(res.pam$silinfo)
## List of 3
## $ widths : num [1:167, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:167] "Finland" "Iceland" "Denmark" "Sweden" ...
## .. ..$ : chr [1:3] "cluster" "neighbor" "sil_width"
## $ clus.avg.widths: num [1:6] 0.5571 0.2517 0.3405 0.0867 0.2638 ...
## $ avg.width : num 0.266
El primero, los widths, es donde esta la información de cada uno de los casos:
# veamos solo algunos.
head(res.pam$silinfo$widths)
## cluster neighbor sil_width
## Finland 1 2 0.6829425
## Iceland 1 2 0.6818799
## Denmark 1 2 0.6763410
## Sweden 1 2 0.6637676
## Australia 1 2 0.6162141
## New Zealand 1 2 0.6122470
Creemos un data frame:
poorPAM=data.frame(res.pam$silinfo$widths)
poorPAM$country=row.names(poorPAM)
Nos interesa sólo sil_width negativos:
poorPAMcases=poorPAM[poorPAM$sil_width<0,'country']
# osea:
poorPAMcases
## [1] "Senegal" "Bhutan" "Sri Lanka"
## [4] "Indonesia" "Bangladesh" "Papua New Guinea"
## [7] "Turkey" "Hong Kong" "Guatemala"
## [10] "Montenegro" "Honduras" "Namibia"
## [13] "Belarus" "Eswatini" "Ethiopia"
La cantidad de paises es:
length(poorPAMcases)
## [1] 15
Podemos hacer lo mismo para las demás estrategias:
# agnes
poorAGNES=data.frame(res.agnes$silinfo$widths)
poorAGNES$country=row.names(poorAGNES)
poorAGNEScases=poorAGNES[poorAGNES$sil_width<0,'country']
#diana:
poorDIANA=data.frame(res.diana$silinfo$widths)
poorDIANA$country=row.names(poorDIANA)
poorDIANAcases=poorDIANA[poorDIANA$sil_width<0,'country']
Ahora que tenemos todos los paises mal asignados, podemos interrogar a los resultados usando teoría de conjuntos, por ejemplo:
intersect(poorAGNEScases,poorPAMcases)
## [1] "Hong Kong" "Bhutan" "Bangladesh" "Honduras" "Guatemala"
setdiff(poorAGNEScases,poorPAMcases)
## [1] "Uruguay" "Mauritius" "South Africa" "Malaysia"
## [5] "Benin" "Georgia" "Ivory Coast" "Fiji"
## [9] "Myanmar" "Niger" "Jordan" "Kuwait"
## [13] "Zimbabwe" "Angola"
setdiff(poorPAMcases,poorAGNEScases)
## [1] "Senegal" "Sri Lanka" "Indonesia"
## [4] "Papua New Guinea" "Turkey" "Montenegro"
## [7] "Namibia" "Belarus" "Eswatini"
## [10] "Ethiopia"
union(poorPAMcases,poorAGNEScases)
## [1] "Senegal" "Bhutan" "Sri Lanka"
## [4] "Indonesia" "Bangladesh" "Papua New Guinea"
## [7] "Turkey" "Hong Kong" "Guatemala"
## [10] "Montenegro" "Honduras" "Namibia"
## [13] "Belarus" "Eswatini" "Ethiopia"
## [16] "Uruguay" "Mauritius" "South Africa"
## [19] "Malaysia" "Benin" "Georgia"
## [22] "Ivory Coast" "Fiji" "Myanmar"
## [25] "Niger" "Jordan" "Kuwait"
## [28] "Zimbabwe" "Angola"
La estrategia basada en densidad sigue una estrategia muy sencilla: juntar a los casos cuya cercanía entre sí los diferencia de otros. Aquí hay un resumen breve del tema:
El algoritmo dbscan requiere dos parametros:
El valor k que se usará es al menos la cantidad de dimensiones. Como se han usado cinco variables, usaremos k=5. Calculemos el epsilon.
Sin embargo, el principal problema es que necesitamos un mapa de posiciones para todos los paises. Eso requiere una técnica que proyecte las dimensiones originales en un plano bidimensional. Para ello usaremos la técnica llamada escalamiento multidimensional:
proyeccion = cmdscale(g.dist,eig=TRUE, k=2,add = T) # k is the number of dim
Habiendo calculado la proyeccción, recuperemos las coordenadas del mapa del mundo basado en nuestras dimensiones nuevas:
# data frame prep:
inputData$dim1 <- proyeccion$points[,1]
inputData$dim2 <- proyeccion$points[,2]
Aquí puedes ver el mapa:
base= ggplot(inputData,aes(x=dim1, y=dim2,label=row.names(inputData)))
base + geom_text(size=2)
Obtengamos los clusters anteriores:
inputData$pam=as.factor(res.pam$clustering)
inputData$agnes=as.factor(res.agnes$cluster)
inputData$diana=as.factor(res.diana$cluster)
Grafiquemos:
# Estimado limites:
min(inputData[,c('dim1','dim2')]); max(inputData[,c('dim1','dim2')])
## [1] -0.6287386
## [1] 0.6412899
limites=c(-0.7,0.7)
base= ggplot(inputData,aes(x=dim1, y=dim2)) + ylim(limites) + xlim(limites) + coord_fixed()
base + geom_point(size=2, aes(color=pam)) + labs(title = "PAM")
base + geom_point(size=2, aes(color=agnes)) + labs(title = "AGNES")
base + geom_point(size=2, aes(color=diana)) + labs(title = "DIANA")
Ahora calculemos usando dbscan:
# euclidea!!
g.dist.cmd = daisy(inputData[,c('dim1','dim2')], metric = 'euclidean')
library(dbscan)
kNNdistplot(g.dist.cmd, k=5)
library(fpc)
db.cmd = dbscan(g.dist.cmd, eps=0.06, MinPts=5,method = 'dist'); db.cmd
## dbscan Pts=167 MinPts=5 eps=0.06
## 0 1 2 3
## border 11 2 0 12
## seed 0 49 33 60
## total 11 51 33 72
Aqui asignamos los clusters a otra columna:
inputData$dbCMD=as.factor(db.cmd$cluster)
Aquí sin texto:
library(ggrepel)
base= ggplot(inputData,aes(x=dim1, y=dim2)) + ylim(limites) + xlim(limites) + coord_fixed()
dbplot= base + geom_point(aes(color=dbCMD))
dbplot
Aquí con mucho texto:
dbplot + geom_text_repel(size=5,aes(label=row.names(inputData)))
Aquí sólo los atípicos:
LABEL=ifelse(inputData$dbCMD==0,row.names(inputData),"")
dbplot + geom_text_repel(aes(label=LABEL),
size=5,
direction = "y", ylim = 0.45,
angle=45,
segment.colour = "grey")
Nota que en esta técnica hay casos que no serán clusterizados.